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Abstract 



We present numerical results on the hydrodynamic stability of coalescing binary 
stars in the first post Newtonian(lPN) approximation of general relativity. We pay 
particular attention to the hydrodynamical instability of corotating binary stars in 
equilibrium states assuming the stiff polytropic equation of state with the adiabatic 
constant T = 3. In previous 1PN numerical studies on corotating binary stars in 
equilibrium states, it was found that along the sequence of binary stars as a function 
of the orbital separation, they have the energy and/or angular momentum minima 
where the secular instability sets in, and that with increase of the 1PN correction, the 
orbital separation at these minima decreases while the angular velocity there increases. 
In this paper, to know the location of the innermost stable circular orbit(ISCO), we 
perform numerical simulations and find where the hydrodynamical instability along 
the corotating sequences of binary sets in. From the numerical results, we found that 
the dynamical stability limit seems to exist near the energy and/or angular momentum 
minima not only in the Newtonian, but also in the 1PN cases. This means that the 
1PN effect of general relativity increases the angular frequency of gravitational waves 
at the ISCO. 



§1. Introduction 

The laser interferometric gravitational wave detectors such as LIGO 0) , VIRGO !\ GEO600 
& and TAMA30oi' ) will be in operation within four years or so. One of the most important 
astrophysical sources of gravitational waves for these detectors is coalescing binary neutron 
stars(BNS's) because gravitational waves emitted in the late inspiraling phase have frequen- 
cies in the sensitive region of these detectors, i.e., from 10Hz to 1000Hz. We will be able to 
know each mass and spin of BNS's if we can obtain an accurate theoretical template for data 
analysis &. For that reason, much theoretical effort has been paid to obtain a theoretical 
template as accurate as possible i\ 

When the orbital separation of BNS's becomes a few times of the neutron star(NS) radius, 
the hydrodynamical effect becomes important. In such a post inspiraling phase, the wave 
form of gravitational waves is expected to be sensitive to the structure of NS such as the 
relation between the radius and the mass of NS. Thus, if gravitational waves from such a 
phase are detected, we may constrain the equation of state(EOS) of NS&&&. In particular, 
it is important to know the innermost stable circular orbit (ISCO), where the binary will 
change the orbit from the stable inspiraling orbit to the merging one. Since its location will 
be sensitive to the structure of NS, it will bring us an important information on the EOS of 
NS. 

The location of the ISCO is determined both by (1) the pure general relativistic(GR) 
effect on two spherical bodies and by (2) the hydrodynamical effect. Recently, Lai, Rasio 
and Shapiro(LRS) pointed out the importance of the latter effect §: They showed that when 
the orbital separation of BNS's becomes sufficiently small, but larger than the approximate 
radius of the ISCO for two spherical bodies ~ 6GM/c 2 where M is the total mass of the 
system, each star of binary is significantly deformed by the tidal force from the companion 
star. Then, an additional tidal field due to the deformation of each star is generated and 
as a result, the circular orbit of BNS's becomes hydrodynamically unstable before two stars 
come into contact even in the Newtonian case if the EOS of NS is stiff enough. 

Such a tidal effect will be sensitive to the structure of NS because the degree of the tidal 
deformation depends on it. Since general relativity plays an important role in determining 
the structure of NS, we can expect that this is also the case in determining the location 
where the hydrodynamical stability due to the tidal force sets in. To know the effect of 
general relativity, as a first step, Shibata0)0) numerically obtained the equilibrium state of 
corotating binary of equal mass in the 1PN approximation and investigated the stability of 
such binaries. The purpose of his works was to clarify the 1PN effect of general relativity 
in BNS's. The main conclusions he obtained are (a) along the equilibrium sequence, there 
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exist the energy and/or angular momentum minima where the secular instability of binary 
sets ini*, and (b) the orbital separation at the minima decreases due to IPN correction, and 
as a result, the orbital angular velocity there increases. 

As LRS showed for the Newtonian case, the secular instability limit does not coincide 
with the dynamical one for corotating binary, and the orbital separation of the latter one 
will locate slightly inside the former one!-*. The ISCO is defined as the dynamical instability 
limit, which is found only by performing numerical simulation. Hence, in order to determine 
the dynamical instability limit along the sequence of corotating binary, we should perform 
numerical simulations and judge the dynamical stability. In this paper we will show that 
the hydrodynamical instability occurs near the energy and/or angular momentum minima 
not only in the Newtonian cases but also in the IPN cases. 

This paper is organized as follows. In §2, we show the IPN hydrodynamic equations to be 
solved which are slightly different from previous ones@: In particular, we change the form 
of the energy equation. In §3, we briefly describe the method for numerical simulation. In 
§4, numerical results are shown: In this paper, we pay attention to binaries of the polytropic 
EOS of the adiabatic index r = 3 for which the energy and angular momentum minima 
for equilibrium sequences were clearly obtained for both the Newtonian and IPN cases in 
previous papers^^. We show numerical results and discuss the dynamical stability of 
binaries around those minima. §5 is devoted to summary. Throughout this paper, G, c, and 
M Q denotes the gravitational constant, the speed of light and the solar mass. 



§2. Basic equations 
In the IPN approximation with the standard gauge condition, the line element is written 



where 
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ds A = g^dx^dx" = —a z c z dt z H — -fiidx l dt + ip 4 "fijdx l dx : ' 



cr 



U 1 fU 2 \ 

a = l -^ + v\~t + X *)+0(c-% (2-1) 

^ = l + § + 0(c~ 4 ), (2-2) 

7tf = * o - + 0(c- 4 ). (2-3) 

Here, (3{ is calculated from 

A = ~\p» + \ WP*i,i + X,i) • (2-4) 
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Equations for gravitational potentials, U#, X*, and x-> are shown below. 
The energy-momentum tensor is 

= p(c 2 + e + ^\ u»u v + Pg» u , 



(2-5) 



where p, e, P and are the baryon density, the internal energy, the pressure, and the four 
velocity, respectively. In the 1PN approximation, components of four velocity are written as 



u° = 1 + I|it; 2 + ^| + 0(CT 4 ), 

l + ^-U,] + 0(c- 
l + \\\v 2 + u\+0(c 



U = - 



u= — 



c 2 2 



v l 1 



Ui = - + ^<p i + v t [ — + 3U*\ +0(c - 5 ), 



c c° 



where 



v % = ^-r and v 2 = v l v % . 



As the EOS, we assume the polytropic one as 

p = (r - l)pe. 

The continuity equation is 

V„(pu") = . 

Using the variable p* = pau°ip 6 , we rewrite the continuity equation as 

dp* d(p*v l ) 



In the 1PN approximation, p* is 



1 /v 



= . 



From the conservation of the energy momentum tensor, 

1 d 



V T 71 — 



aip 6 [pc 2 + pe + P I v^Vby 



1 



(2-6) 
(2-7) 

(2-8) 
(2-9) 

(2T0) 

(2-H) 
(2-12) 

(2-13) 
(2-14) 



0, 



(2-15) 
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we obtain eqautions for the Euler and the entropy conservation as 

H — — = -aip P ti - p*cm u a ; + p*Uj(3 J ^ - -^p u j u kT ,i 



dt dxi 
<9e* d(e*v j ) 

dt dx j 



0. 



where 



Wi= 4 1 + 





p 








p 




p 








p 



>u 



e*= (pe) 1/r au ^. 
In the 1PN approximation, the Euler equation becomes 



dt 



1 + 



2U* 



X 



\Pa + P* 



+ 0(c- 4 ) , 



1 



3v 2 



P 



where 



£ + 



P v 2 



3U.\ } + % + 0(c- 4 ) 



The pressure is calculated from 



P=(r-l)ef 1 + 



1 /V 



(2-16) 
(2-17) 

(2-18) 

(2-19) 
(2-20) 



! + —[£ + - + — -U* \ W-t 



(2-21) 



(2-22) 



(2-23) 



We notice that the form of the 1PN hydrodynamic equations, in particular the entropy 
equation, shown above are different from those of Blanchet, Damour, and Schafer0) although 
in both formalisms, every 1PN term is taken into account consistently. Thus, equations 
adopted in this paper are different from those used in a previous paper©. 

Equations for various potentials [/*, X*, P*i and x are derived from the Einstein equation 

as 



AU* = —AnGp*, 
AP„ = -AirGp*v\ 



AX, = AiiGp* 
Ax = 4nGp*v l x % 



3P 

— - U, + £ + — 

2 p 



(2-24) 
(2-25) 

(2-26) 

(2-27) 
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Note that the definitions of the gravitational potentials are based on p # , not on p. Thus, 
they are different quantities from those defined in obtaining the equilibrium state 0)0) in 
which we need to use the gravitational potentials defined with respect to p. (For example, 
the Newtonian potential U which satisfies AU = —Airp differs from [/*. ) 

We use the following forms of the continuity, Euler, and the entropy equations as the 
basic equations in the Newtonian case: 

d(pv l ) | (9(pt>V) _ dP dU 
dt dxi dx % dx l ' 

where e = (pe) 1 ^ = [P/(r — l)] 1 ^. Thus, we change the form of the energy equation from 
that adopted in previous papers0)0). We emphasize that this choice improves accuracy on 
the conservation of the energy and angular momentum in numerical simulation compared 
with our previous works. 

Before closing this section, we define the following quantities correct up to 1PN order: 

• Conserved mass and Newtonian mass 



M ± 



J pJV and M = J pdV. (2-31) 



Note that for the Newtonian case, M* = M. 
The energy: In the definition based on p, 

+2eU- ^U 2 ^j \dV, (2-32) 
while in the definition based on p*, 

-eU m + ±U^dV 
£*. (2-33) 
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• The angular momentum around the rotation axis : In the definition based on p, 



J 




{-yv x + xv y ) ll + \\v 2 + W + e+-\\ + \{-y(3 x + x(3 y ) 



dV, (2-34) 



while in the definition based on p*, 

J = J p*(—yu x + xu y )dV = J*. (2-35) 

Note that we adopt the z-axis as the rotation axis in this paper. Here, we also notice 
that numerical values E and E* or J and J* are, respectively, slightly and system- 
atically different because some terms of 0(c~ 4 ) are implicitly included in E* and J*. 
However, the difference is small and not important in the following discussions. In the 
following, we use J* as the angular momentum in numerical simulations, while we use 
the definition based on p as the angular momentum of the equilibrium sequence. 
Center of mass of each star 

IVlif J each star 

From Xg, we define the coordinate separation of the orbital radius as r g = 2(x l g x g ) l l 2 . 
The quadrupole moment of system 

I.. = J prfjPdV. (2-37) 

We also define 

Irr = Ixx + Iyy, and -t^ = 1^ — -5^1^, (2-38) 
where 5^ is the Kronecker's delta. 



§3. Numerical method 

In numerical computation, we solve (1) the 1PN hydrodynamic equations and (2) the 
Poisson equations for various potentials. In the following, we mention the numerical methods 
we adopt briefly. 

Our treatment for the advection term in the hydrodynamic equations is the same as that 

adopted in previous papers 0)0)0). Although the method is the same, the energy equation 

we adopt is different from the previous one. As a result, accuracy on the conservation of 

the angular momentum and the energy of the system is significantly improved.0 As for the 

*' We here note that using the entropy equation adopted above, the entropy is conserved. If the shock 
is formed, however, the entropy should not be conserved. Hence, we must introduce the artificial viscosity 
term in the Euler and entropy equations to express the shock accurately. In this paper, however, the shock 
formation does not play an important role, so that we do not add such terms. 
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time evolution, we use the second order Runge-Kutta method. Thus, the numerical code 
for solving the hydrodynamic equations is second order accurate both in time and space. 
(Actually, we checked, for example, that accuracy of the conservation of the total energy is 
second order convergent changing the grid number.) 

Since we use the stiff EOS, density near the surface of each star decreases steeply. With 
a finite resolution of numerical grid, one cannot accurately represent the density decline 
near the surface. As a result, the pressure gradient often happens to become too steep, 
and the velocity near the surface becomes unphysically too large. To avoid that, when the 
density p* is less than p crit ~2x 10~ 4 p max , where p max is the maximum of p*, we artificially 
suppress the linear momentum (or v l for the Newtonian case) near the surface by a factor 
p*j (p* + /Pmax); where / is a small factor and we set as ~ 4 x 10~ 5 . This technique has been 
used in our series of papers 0)0)©. 

To solve 6 Poisson equations for U*, X*, P* x , P* y , P* 2 , and x-> we use the ICCG method 
as we have done in previous papers© 0)0) imposing the boundary conditions at outer grids 
as 



GM, 




P 

P*x ^ / p*xv x dV + ^ / p*yv x dV + 0(r~ 4 ), (3-3) 
r J v J 

P *y^^r[ p* xyVdv + ~r / p*v vVdv + c(r~ 4 ), (34) 

J p*zv z dV + 0(r- 4 ), (3-5) 



(3-1) 

U*+e + - )dV + 0(r~ 3 ), (3-2) 



Gz 



/ + 0(r- 3 ). (3-6) 

r J 

Thus, we use three types of numerical implementations for solving Poisson equations for (Z7#, 
X*, x), (P*xi P*y)i an d P*zi respectively. Accuracy of these Poisson solvers has been checked 
to be less than 0.1% and a result for the case of a small grid number is shown in0). 

Throughout this paper, we assume the symmetry with respect to the equatorial plane. 
Thus, we also impose boundary conditions at the equatorial plane as 

p*{—z) = p*{z), e*(-z) = e*(z), v A {-z) = v A (z), 

v *(- z ) = -v z {z), U*{-z) = U*{z), X,(-z) = X„(z), 

X (-z) = X (z), P*a(~z) = P* A (z), P mz (-z) = -P* z (z), (3-7) 

where A = x or y. Grid number we adopt is (N x , N y , N z ) = (141, 141, 71), and grid spacing 
is chosen in order that the grid covers the minor axis of each star at initial state by at 



8 



least 15 grid points. Numerical computations were performed on FACOM VX4 in the data 
processing center of National Astronomical Observatory in Japan. 

FACOM VX4 has four parallel processors, so that four different procedures can be done 
on separate processors at the same time. To solve the Poisson equations in the 1PN case, we 
make use of this property: In each time step, we solve the Poisson equations for U* and X* 
on processor 1, for P* x and P* y on processor 2, for P* z on processor 3, and for x on processor 
4, while the hydrodynamic equations are solved only on processor 1. Since computation of 
solving the Poisson equations is the most time-consuming part in the numerical implemen- 
tations, the CPU time required reduces by a factor of ~ 3 in the 1PN case using this simple 
technique. On the other hand, we only use one processor in the Newtonian case because 
only one Poisson equation is required to be solved. In a typical computation for one model, 
20000 ~ 30000 time steps were necessary, and computation of 20000 time steps took about 
30 CPU hours for the Newtonian case, and about 65 CPU hours for the 1PN case. 



§4. Numerical results 

We assume r = 3 as the adiabatic index because we are interested in BNS's which are 
composed of a stiff EOS, i.e., r = 2 ~ 30- 1 . As initial conditions of binary, we adopt 
equilibrium states of binary corotating around the z-axis. When we obtain such a solution, 
we set the EOS as 

P = Kp r , (4-1) 

where K is the polytropic constant. Following a previous paper 0, for convenience, we set 
K as 

K = 2.52AGr 5 M-\ (4-2) 

Then, the radius of each star of binary in the Newtonian limit is r o (M/M ) a2 = ao- In the 
following, M* denotes the conserved mass of each star of binary and is set as 1.4M Q . In this 
paper, we consider two cases, the Newtonian case(i.e., c — > oo in the PN approximation) 
and the 1PN case of r = 40R S where R s = GM Q /c 2 . We note that for r = 40i? s , the 
compactness of each star is C s = GM^/a^c 2 ~ 0.033, where a* = r o (M„/M ) o ' 2 (note that 
a* = ao in the Newtonian case), and this value is too small to represent realistic BNS's 
which we are most interested in. However, the PN approximation is valid only for weak 
gravity, i.e., for small compactness, and it is meaningless to apply it for large compactness. 
We emphasize that the main aim of this paper is to solidly know the 1PN effect which will 
play an important role in BNS's. 
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In figs.l, we show the energy and the angular momentum of corotating binary in equilib- 
rium states as a function of the orbital separation r g /a* for the 1PN case of tq = 40R S (filled 
circles) and the Newtonian case(open circles). As shown in previous papers^ 0), there exist 
the minima along the sequences of the energy and angular momentum, and the separation 
at the minima becomes small with increase of the 1PN effect: For the Newtonian case, the 
separation is r g ~ 3.08a , and for the 1PN case of r = 40-R s , it is r g ~ 2.94a*. 

We will perform numerical simulations using some of such equilibrium states near the 
energy and/or angular momentum minima as initial conditions. Among equilibrium states 
along each sequence, we choose four states as initial conditions: they are N5, N6, N7, and 
N8 along the Newtonian sequence, and PN5, PN6, PN7, and PN8 along the 1PN sequence, 
respectively. Models N6 and PN6 are equilibrium states at the nearest location from the 
energy and/or angular momentum minima. 

In figs. 2 (a) and (b), we show the evolution of Irr for the Newtonian and 1PN cases, 
respectively. Hereafter, we use the unit of time as 4n^Jal/GM*, and in this unit, the orbital 
periods of equilibrium states are ~ 1.82 for N5, 1.88 for N6, 1.95 for N7, 2.03 for N8, 1.74 for 
PN5, 1.80 for PN6, 1.87 for PN7, 1.95 for PN8, respectively. On the other hand, the vertical 
axis is shown in units of M*a*(i.e., / = Irr/ '(M*al)). This means that when I < 2, the 
binary begins to merge. Hereafter, solid, dotted, dashed and dotted-dashed lines are used 
to show quantities for models N5(PN5), N6(PN6), N7(PN7), and N8(PN8), respectively. 

Figs. 2 show that equilibrium models N5, N6, PN5 and PN6 seem hydrodynamically 
unstable; two stars merge within about one orbital period. On the other hand, models N7, 
N8, PN7, and PN8 seem hydrodynamically stable because in these binaries, two stars do 
not come into merging in the dynamical time scale(i.e., the orbital period). We note that 
irrespective of the hydrodynamical stability, the orbit of binaries shrinks due to dissipation of 
the angular momentum. This occurs because of a finite resolution(i.e., finite accuracy) both 
in time and space in the numerical simulation, and is the reason why I decreases gradually 
for stable binaries such as N7, N8, PN7, and PN8. 

In fig. 3, we show the angular momentum as a function of time for the Newtonian and 
1PN cases( note that for the 1PN case, the angular momentum is J*). Here, the unit 
of the angular momentum is G 1 l 2 M^ 2 a 1 J 2 , and the unit of the time is the same as that 
in figs. 2. Notice that a small fraction of matter carrying the large angular momentum is 
ejected outside the numerical grids in model PN5 for t ^ 2.6, and in model N5 for t ^ 2.9 so 
that the angular momentum is decreasing. Fig. 3 shows that accuracy on the conservation of 
the angular momentum is fairly good in all simulations, but it is not perfect: In one orbital 
period, the angular momentum is lost by about ~ 0.1%. 

Even for the angular momentum loss, equilibrium models N7, N8, PN7, and PN8 do 
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not proceed to merging in the dynamical time scale. Thus, these binaries seem obviously 
stable. On the other hand, we cannot completely deny the possibility that a small loss of 
the angular momentum by such numerical effects leads models N5, N6, PN5 and PN6 to 
unstable orbits. However, the loss of the angular momentum is not so large up to the contact 
of these binaries. According to an approximate analysis by LRS(see, for example, fig. 1 in&), 
the hydrodynamically stable corotating binary will move to another stable binary which 
is not corotating state even when it slightly loses the angular momentum (unless its loss is 
so large). Hence, we judge that the hydrodynamical instability will really occur when the 
orbital separation of binary becomes as small as that of N5, N6, PN5 and PN6 or slightly 
smaller than these, i.e. very near or slightly inside the energy and/or angular momentum 
minima along the sequences of corotating binaries, r g ~ 3.08a for the Newtonian case, and 
r g ~ 2.94a, for the 1PN case of C s ~ 0.033. 

We note the following fact: As for the Newtonian study on the determination of the 
dynamically stability limit for the r = 3 case, there have been two independent works by 
Rasio and Shapiro(RS)§) and New and Tohline(NT)@), and their results slightly disagree: 
According to RS, the dynamical instability along the corotating sequence sets in at r g ~ 
2.97ao, while according to NT, it sets in at r g ~ 3.08ao- The hydrodynamic code which RS 
used in their numerical simulations is a SPH one, while that used by NT is the Eulerian 
one. In our present simulation, we use the Eulerian code, and this may be the reason why 
the present results support the work by NT. Currently, however, the reason for this slight 
difference is not clear. 

In fig. 4, we show the maximum density of the system as a function of time for the 
Newtonian and 1PN cases, respectively. The density is shown in units of M ! „(47ra2/3) _1 . 
Fig.4 shows that apart from a small oscillation, the maximum density is almost constant for 
the case where two stars of the binary do not merge. We notice that the maximum density 
for the 1PN configuration are larger than that for the Newtonian one for a given angular 
frequency. This is a manifestation of the GR effect by which each star of binary is forced 
to be more compact than that in the Newtonian case. It should be also noted that along 
not only the Newtonian equilibrium sequence, but also the 1PN one, the central density 
decreases with decrease of the orbital separation©©. 

In figs. 5 and 6, we show the evolution of the density profile in the equatorial plane for 
models N7 and PN7. Note that x- and y-axes are shown in units of a*, outermost grid 
is placed at 4.09 in this unit, and the orbital rotation is counterclockwise. Figs. 5 and 6 
show that for more than one orbital period, the density profile does not change so much 
apart from a small deformation in the outer region of each star. Thus, in our present 
numerical simulation, the density profile of equilibrium state can be kept fairly well for 
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hydrodynamically stable models. 

In figs. 7 and 8, we show the evolution of the density profile in the equatorial plane for 
models N5 and PN5, i.e., dynamically unstable models, respectively. In both Newtonian 
and 1PN cases, the evolution is similar to that obtained by RS©: In about one orbital 
period, two stars merge due to the hydrodynamic instability. After they merge, spiral arms 
are formed and they transport the angular momentum outward. In the final stage, the 
spiral arms widen and merge together. In the central region, an ellipsoidal core is formed. 
Since radiation reaction of gravitational waves is not taken into account in this paper, the 
ellipsoidal figure will be kept forever. 

Finally, we show the wave form and the luminosity of gravitational waves for merging 
phase of binary. Although we are able to calculate the 1PN wave form and luminosity 
using the Blanchet, Damour, and Schafer formalism0\ in this paper we simply show the 
Newtonian ones. In figs. 9(a) and (b), we show + and x modes of gravitational waves 
observed along the z-axis for models PN5 (solid lines) and N5 (dotted lines). + and x modes 
are defined as 

G d 2 ( T T \ , , 2G d 2 r 

r ° h + = ZaZu2 I l ™ ~ !yy I ' and r ° h * = —^i 1 ^ ( 4 ' 3 ) 



where r denotes the distance from the source to an observer on the z-axis. In fig. 9(c), we 
show the luminosity defined as 

dE G d?-t;j d 3 ^-^ . , ,. 

— = -. (44) 

dt 5c 5 dt 3 dt 3 1 ; 

Since the orbital period of model PN5 is shorter than that of N5, the wave length of 
gravitational waves before merging for PN5 is shorter than that of N5. Also, the each star of 
model PN5 is more compact than that of N5, so that the maximum amplitude of the lumi- 
nosity for PN5 is larger than that for N5. Although these differences exist between PN and 
Newtonian results, overall shapes of the wave form and luminosity in the PN approximation 
are essentially the same as those in the Newtonian case (see also figs. 10 and 11 in ref . Elt ) : As 
the orbit shrinks due to the hydrodynamical instability, the amplitude and the luminosity 
of gravitational waves gradually rise. When two stars of the binary come into contact, they 
become the maximum, and after that, they gradually fade away. Since the merged object 
becomes an ellipsoid in the final state and radiation reaction of gravitational waves is not 
included in this simulation, the wave form will settle down to a stationary one. 

Before closing the section, we notice that evolution of the merged object and gravitational 
waves from it considerably depend on the velocity field of pre-merging state of binary: In the 
case where each star of binary does not have spin such as the irrotational Roche-Riemann 
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binary which is regarded as one of the most realistic situations for BNSE-P, we have learned 
from numerical simulations that the large spiral arm is not formed Also, when radi- 

ation reaction of gravitational waves is included and as a result, pre-merging binary has an 
approaching velocity, the merged object radially oscillates and second and third peaks of the 
gravitational wave luminosity, which in the present result do not appear, can be seen^H*©. 
Thus, we should keep in mind that the present results on the evolution of merged object 
and gravitational waves from it are inherent for evolution of dynamically unstable corotating 
binary without radiation reaction. 



§5. Summary 

In this paper, we have shown numerical results of the hydrodynamical stability of coro- 
tating binary in equilibrium states in the 1PN approximation of general relativity. We found 
that the hydrodynamical instability sets in near the energy and/or angular momentum min- 
ima along the equilibrium sequence. This means that we may approximately determine the 
location of the ISCO from the energy and/or angular momentum minima along sequences 
for corotating binary. 

As shown in previous papers, the angular velocity at those minima changes with the first 
GR correction approximately as 

^ = ^n(i + Cpn^), (5-1) 

where i?N is the angular velocity at the minima for the Newtonian case and Cpn is a constant 
- 1.1 for r = 3 and ~ 2.5 for F = 1. For the polytropic EOS of F = 30 } and 20, 



(1 N ~ (0.267 ±0.002) J forT = 3, (5-2) 



~ (0.297 ± 0.002V— ^ for r = 2, (5-3) 
V a o 

and the frequency of gravitational waves is calculated as 

3/2 

h = —^ 6301-^—1 ( — I Hz forr = 3, (5-4) 




M \ 


1/2 


( 15km 


IAMqJ 


1 1 


^ ao 


M \ 


1/2 


I 15km 


1.4M J 


1 1 


1 a 



3/2 

700 1 7T7 I ! ~~- L I Hz for F = 2 ' ( 5 ' 5 ) 
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where M and ao denote the Newtonian mass and radius of each star of binary, respectively. 
Since the EOS of NS has the adiabatic index of 2 ~ 3, /n may be approximated by 

/ \ 1/2 / \ 3/2 

Thus, in the case of corotating binary, the frequency of gravitational waves at the ISCO for 
a realistic BNS of radius 10~15km(i.e., C s = 0.14 ~ 0.2) will become about 150 ~ 300Hz 
higher due to the 1PN correction of general relativity. 

We performed numerical simulations only including the 1PN correction, i.e., the lowest 
order GR correction, in this paper, but it is not sufficient in order to obtain the ISCO for 
realistic BNS's accurately because BNS's are highly GR objects(C s ~ 0.2). We need to 
take into account the fully GR term to obtain an accurate result. We, however, emphasize 
that the conclusion in this paper will hold for binaries of small compactness(say C s ^ 0.05). 
This means that when we perform fully GR simulation, the numerical code must reproduce 
the present result. Thus, the present paper will be a useful guideline in checking fully GR 
calculations. 
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Figure captions 

Fig.l The energy(a) and the angular momentum(b) of corotating binary in equilibrium states 
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as a function of the orbital separation r 9 /a* for the 1PN case of r$ = 40-R s (filled circles) 
and the Newtonian case(open circles). The energy and the angular momentum are shown 
in units of GM 2 /a* and G x ^ 2 mV 2 aj 2 , respectively. 

Fig.2 Normalized moment Irr/(M*o,1) as a function of time for the Newtonian case(a) and 
the 1PN case of C s ~ 0.033(b). Units of time is 4nJal/GM,. Solid, dotted, dashed 
and dotted-dashed lines are for models N5(PN5), N6(PN6), N7(PN7), and N8(PN8), 
respectively. 

Fig.3 The angular momentum as a function of time. Units of the angular momentum and the 
time are G 1 ! 2 M^ 2 aV 2 and Air^Jal/GM*, respectively. Solid, dotted, dashed and dotted- 
dashed lines are for models N5(PN5), N6(PN6), N7(PN7), and N8(PN8), respectively. 

Fig.4 The same as fig.3, but for the maximum density. Units of the density is M^Aita^/SY 1 . 

Fig.5 Time evolution of the density(p) profile in the equatorial plane for model N7. x and 
y-axes are shown in units of a*(= a ). Note that outermost grid is placed at 4.09 in 
this units, and the orbital rotation is counterclockwise. Density contours (solid lines) are 
spaced with intervals of p max ,o/10, where p max ,o is the maximum density at the initial state. 
Dotted lines denote contour of p ma x,o/100. 

Fig. 6 The same as fig.5, but for of model PN7. 

Fig. 7 The same as fig.5, but for model N5. 

Fig. 8 The same as fig.5, but for of model PN5. 

Fig. 9 Wave forms r Q /i + (a) and r h x (b), and the luminosity(c) of gravitational waves for 
models PN5(solid lines) and N5(dotted lines). Units of the amplitude of r Q h + and r a h x , 
the luminosity, and time are G 2 M 2 /c 4 a*, c 5 /G(GM*/a*c 2 ) 5 and AirJal/GM^, respectively. 
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